home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Shareware Grab Bag
/
Shareware Grab Bag.iso
/
090
/
byt1186b.arc
/
RUNGKUT.LBR
/
RKCONST.FOR
< prev
next >
Wrap
Text File
|
1986-04-11
|
6KB
|
237 lines
$NOFLOATCALLS
$NODEBUG
$STORAGE:2
c*************************************************************************
subroutine const(ist,imeth,idef)
c**** subroutine supplies constants used in the runge-kutta
c**** method,in the format required by RKP.
implicit double precision (a-h,p)
dimension al(12),b(12,12)
common /coeffs/al,b,a1,a3,a4,a5,a6,a7,a8,a9,a10,a11,
& a12,a13,er1,er3,er4,er5,er6,er7,er8,
& er9,er10,er11,er12,er13,inum
common /cons/hmin,hmaxl,hfact,hlim,power
do 10 j=1,12
do 5 i=1,12
b(i,j)=0.d0
5 continue
al(j)=0.d0
10 continue
if(imeth.EQ.1) then
c**** Fourth order runge-kutta method specified by
c**** Felberg, E., COMPUTING, 6(1970)p61-71
c write(*,*) 'Fourth order Runge-Kutta-Fehlberg method'
al(1)=1.d0/4.d0
al(2)=3.d0/8.d0
al(3)=12.d0/13.d0
al(4)=1.d0
al(5)=1.d0/2.d0
b(1,1)=1.d0/4.d0
b(1,2)=3.d0/32.d0
b(2,2)=9.d0/32.d0
b(1,3)=1932.d0/2197.d0
b(2,3)=-7200.d0/2197.d0
b(3,3)=7296.d0/2197.d0
b(1,4)=439.d0/216.d0
b(2,4)=-8.d0
b(3,4)=3680.d0/513.d0
b(4,4)=-845.d0/4104.d0
b(1,5)=-8.d0/27.d0
b(2,5)=2.d0
b(3,5)=-3544.d0/2565.d0
b(4,5)=1859.d0/4104.d0
b(5,5)=-11.d0/40.d0
er1=1.d0/360.d0
er3=-384.d0/12825.d0
er4=-41743.d0/1429560.d0
er5=1.d0/50.d0
er6=2.d0/55.d0
a1=16.d0/135.d0
a2=0.d0
a3=6656.d0/12825.d0
a4=28561.d0/56430.d0
a5=-9.d0/50.d0
a6=2.d0/55.d0
inum=6
ist=4
if(idef.eq.0) hmaxl=0.5d0
hlim1=3.d0
end if
if(imeth.EQ.2) then
c**** Fifth order runge-kutta method specified by
c**** Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 5)
c write(*,*) 'Fifth order Runge-Kutta-Verner method'
al(1)=1.d0/18.d0
al(2)=1.d0/6.d0
al(3)=2.d0/9.d0
al(4)=2.d0/3.d0
al(5)=1.d0
al(6)=8.d0/9.d0
al(7)=1.d0
b(1,1)=1.d0/18.d0
b(1,2)=-1.d0/12.d0
b(2,2)=1.d0/4.d0
b(1,3)=-2.d0/81.d0
b(2,3)=4.d0/27.d0
b(3,3)=8.d0/81.d0
b(1,4)=40.d0/33.d0
b(2,4)=-4.d0/11.d0
b(3,4)=-56.d0/11.d0
b(4,4)=54.d0/11.d0
b(1,5)=-369.d0/73.d0
b(2,5)=72.d0/73.d0
b(3,5)=5380.d0/219.d0
b(4,5)=-12285.d0/584.d0
b(5,5)=2695.d0/1752.d0
b(1,6)=-8716.d0/891.d0
b(2,6)=656.d0/297.d0
b(3,6)=39520.d0/891.d0
b(4,6)=-416.d0/11.d0
b(5,6)=52.d0/27.d0
b(1,7)=3015.d0/256.d0
b(2,7)=-9.d0/4.d0
b(3,7)=-4219.d0/78.d0
b(4,7)=5985.d0/128.d0
b(5,7)=-539.d0/384.d0
b(7,7)=693.d0/3328.d0
er1=33.d0/640.d0
er3=-132.d0/325.d0
er4=891.d0/2240.d0
er5=-33.d0/320.d0
er6=-73.d0/700.d0
er7=891.d0/8320.d0
er8=2.d0/35.d0
a1=57.d0/640.d0
a3=-16.d0/65.d0
a4=1377.d0/2240.d0
a5=121.d0/320.d0
a7=891.d0/8320.d0
a8=2.d0/35.d0
inum=8
ist=5
if(idef.eq.0) hmaxl=1.0d0
hlim1=4.d0
end if
if(imeth.EQ.3) then
c**** Seventh order runge-kutta method specified by
c**** Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 7)
c write(*,*) 'Seventh order Runge-Kutta-Verner method'
al(1)=1.d0/4.d0
al(2)=1.d0/12.d0
al(3)=1.d0/8.d0
al(4)=2.d0/5.d0
al(5)=1.d0/2.d0
al(6)=6.d0/7.d0
al(7)=1.d0/7.d0
al(8)=2.d0/3.d0
al(9)=2.d0/7.d0
al(10)=1.d0
al(11)=1.d0/3.d0
al(12)=1.d0
b(1,1)=1.d0/4.d0
b(1,2)=5.d0/72.d0
b(2,2)=1.d0/72.d0
b(1,3)=1.d0/32.d0
b(3,3)=3.d0/32.d0
b(1,4)=106.d0/125.d0
b(3,4)=-408.d0/125.d0
b(4,4)=352.d0/125.d0
b(1,5)=1.d0/48.d0
b(4,5)=8.d0/33.d0
b(5,5)=125.d0/528.d0
b(1,6)=-1263.d0/2401.d0
b(4,6)=39936d0/26411d0
b(5,6)=-64125.d0/26411d0
b(6,6)=5520.d0/2401.d0
b(1,7)=37.d0/392.d0
b(5,7)=1625.d0/9408.d0
b(6,7)=-2.d0/15.d0
b(7,7)=61.d0/6720.d0
b(1,8)=17176.d0/25515.d0
b(4,8)=-47104.d0/25515.d0
b(5,8)=1325.d0/504.d0
b(6,8)=-41792.d0/25515.d0
b(7,8)=20237.d0/145800.d0
b(8,8)=4312.d0/6075.d0
b(1,9)=-23834.d0/180075.d0
b(4,9)=-77824.d0/1980825.d0
b(5,9)=-636635.d0/633864.d0
b(6,9)=254048.d0/300125.d0
b(7,9)=-183.d0/7000.d0
b(8,9)=8.d0/11.d0
b(9,9)=-324.d0/3773.d0
b(1,10)=12733.d0/7600.d0
b(4,10)=-20032.d0/5225.d0
b(5,10)=456485.d0/80256.d0
b(6,10)=-42599.d0/7125.d0
b(7,10)=339227.d0/912000d0
b(8,10)=-1029.d0/4180.d0
b(9,10)=1701.d0/1408.d0
b(10,10)=5145.d0/2432.d0
b(1,11)=-27061.d0/204120.d0
b(4,11)=40448.d0/280665.d0
b(5,11)=-1353775.d0/1197504.d0
b(6,11)=17662.d0/25515.d0
b(7,11)=-71687.d0/1166400d0
b(8,11)=98.d0/225.d0
b(9,11)=1.d0/16.d0
b(10,11)=3773.d0/11664.d0
b(1,12)=11203.d0/8680.d0
b(4,12)=-38144.d0/11935.d0
b(5,12)=2354425.d0/458304.d0
b(6,12)=-84046.d0/16275.d0
b(7,12)=673309.d0/1636800.d0
b(8,12)=4704.d0/8525.d0
b(9,12)=9477.d0/10912.d0
b(10,12)=-1029.d0/992.d0
b(12,12)=729.d0/341.d0
er1=-1.d0/480.d0
er6=-16.d0/375.d0
er7=-2401.d0/528000.d0
er8=2401.d0/132000.d0
er9=243.d0/14080.d0
er10=-2401.d0/19200.d0
er11=-19.d0/450.d0
er12=243.d0/1760.d0
er13=31.d0/720.d0
a1=31.d0/720.d0
a6=16.d0/75.d0
a7=16807.d0/79200.d0
a8=16807.d0/79200.d0
a9=243.d0/1760.d0
a12=243.d0/1760.d0
a13=31.d0/720.d0
inum=13
ist=7
if(idef.eq.0) hmaxl=2.d0
hlim1=5.d0
end if
power=1.d0/dble(ist+1)
hfact=0.5d0**(1.d0/dble(ist-1))
hlim=hlim1/hfact
if(idef.eq.0) hmin=1.d-06
c write(*,*)' Constants evaluated'
return
end